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Abstract 

We investigate the application of the discontinuous Petrov-Galerkin (DPG) finite el- 
ement framework to stationary convection-diffusion problems. In particular, we demon- 
strate how the quasi-optimal test space norm can be utilized to improve the robustness 
of the DPG method with respect to vanishing diffusion. We numerically compare coarse- 
mesh accuracy of the approximation when using the quasi-optimal norm, the standard 
norm, and the weighted norm. Our results show that the quasi-optimal norm leads to 
more accurate results on three benchmark problems in two spatial dimensions. We ad- 
dress the problems associated to the resolution of the optimal test functions with respect 
to the quasi-optimal norm by studying their convergence numerically. In order to facilitate 
understanding of the method, we also include a detailed explanation of the methodology 
from the algorithmic point of view. 

Keywords: convection-diffusion; finite element method; discontinuous Petrov-Galerkin; 
numerical stability; robustness; automatic stabilization technique 



1 Introduction 

The purpose of this paper is to study the application of the discontinuous Petrov-Galerkin finite 
element methodology to convection-diffusion problems of the form 

— eAu + a ■ V« = f in f2 

u = g on oil 

where Q is a bounded two-dimensional domain, a is the flow velocity, e is the diffusion parameter 
and / represents a source term. The function g determines the values of the solution u on the 
boundary of the domain denoted by dfl. This equation is a basic model for various transport 
processes in nature and engineering. Typical applications include the movement of substances 
(dissolved nutrients, toxins) in the environment. Problem (fTl) may also be viewed as a stepping 
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stone for designing numerical methods for more complex problems such as the Navier-Stokes 
equations for compressible and incompressible flows. 

Construction of numerical solution schemes for that can approximate convection dom- 
inated flows is known to be an inherently difficult task. Many scientists and engineers have 
addressed the problem and designed algorithms based on finite difference, finite volume, finite 
element and spectral approaches. In the context of the finite element method, the standard 
Ritz-Bubnov-Galerkin formulation is known to be deficient and various alternative formulations 
have been proposed in the literature, see for instance jZHIDl El EH E2 EZrES E2 Ell ESI EH] ■ 

The success of the traditional Galerkin finite element method in most structural problems 
is based on its best approximation property. This means that the difference between the 
finite element solution and the exact solution is minimized with respect to a norm induced by 
the strain energy product. Numerical issues arise when the best approximation property is 
lost. This happens, for instance, when the standard Galerkin method is applied to convective 
transport problems. In these problems, the system matrix associated to convection is not 
symmetric and numerical solutions tend to show spurious, non-physical oscillations unless the 
finite element mesh is heavily refined. It should be noted that the accuracy of the finite element 
method may also degrade in problems with symmetric system matrices. Thin-body problems of 
solid mechanics and wave propagation problems constitute typical examples. In these problems, 
the standard Galerkin method leads to non-physical solutions suffering from the phenomena 
of numerical locking (small thickness) and numerical dispersion (high wave number) unless an 
over-refined mesh is used. 

To save computer resources, alternative formulations which produce reasonable approxima- 
tions regardless of the mesh size have been developed. In the context of convective transport 
problems, a well-known method is the Streamline Upwind Petrov-Galerkin (SUPG) formulation, 
see [HI EE] • The method belongs to a larger class of stabilized finite element methods reviewed, 
for example, in [T2J ES]- The stabilized finite element methods are geared mainly towards the 
classical /i-version of the finite element method although some higher order formulations have 
been studied in pfl UU El EH] • 

A new, general finite element framework has been developed by Demkowicz and Gopalakr- 
ishnan in [T5HTT] . This variational framework is of discontinuous Petrov-Galerkin (DPG) type 
and can be utilized, for example, to symmetrize non-symmetric problems by using the concept 
of optimal test functions similar to the early methodology developed by Barrett and Morton 
in [3], see also [211 E2 EH1 EE, ES]. More precisely, if a fully discontinuous finite element space 
is used for the test functions, as in the DPG method of Bottasso et al. [5], the optimal test 
functions can be approximated locally in an enriched finite element space. The same procedure 
can also be used to compute local a posteriori error estimates to guide automatic adaptive mesh 
refinement as demonstrated in [T^] in the context of convection-dominated diffusion problems 
and in [UJ in the context of shell boundary layers. 

Previous work in the DPG framework applied to Problem Q, has utilized adaptive mesh 
refinements [TH [19] to fully resolve fine scales (boundary and interior layers) in the solution. 
However, in engineering applications this is not always possible. In this work, we focus on the 
coarse-mesh accuracy of the DPG method. That is, we investigate how the effect of fine scales 
to the solution can be taken into account without resolving them by the trial space mesh. It 
should be noted that coarse-mesh accuracy is valuable even in the context of adaptive mesh 
refinement to avoid unnecessary mesh refinement. 

We complement the earlier works [TBI ITU] by resolving the stable test function space with 
respect to the optimal (quasi-optimal) test space norm. The corresponding DPG method is 
expected to be uniformly stable with respect to vanishing diffusion, but the local problems 
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for the test functions become singularly perturbed and therefore more difficult to solve. The 
main contribution of our study is to show how this problem can be approached by utilizing a 
carefully designed element sub-grid discretization. The present work is a continuation of our 
recent work [I2J 03] and extends the earlier analysis in one space dimension to two dimensions. 

The paper is structured to be self-contained and focused on bridging the theoretical frame- 
work to a practical understanding. Section [2] lays down the abstract functional analytic setup 
of the method. Then in Section [3j we specialize the theoretical framework to Problem and 
introduce different test space norms. In Section [4], we re-explain the DPG framework from a 
purely algorithmic point of view in order to highlight the practical implications of the different 
choices of the test space norm. We provide details on how to construct the linear algebraic sys- 
tems constituting the local optimal test functions. Numerical results are presented in Section 
[5] and the paper ends with conclusions and remarks in Section |6j 

2 Automatically Stable DPG Framework 
2.1 Background 

The starting point of the formulation is an abstract variational statement of the form: Find 
uGW such that 



where the test function space W and the trial function space U are assumed to be Hilbert 
spaces under inner products (•, -)w and (•, -)u with the corresponding norms || • ||yv an d || • \ \u- 
We shall assume that the spaces are real because our application, Problem Q, deals only with 
real- valued functions. 

If the spaces W and U are the same and B induces an inner product, then the standard 
Galerkin method delivers the best approximation of the exact solution in the corresponding 
norm known as the energy norm. In the case of Problem ([I]) the conventional weak formulation 
is associated with the bilinear and linear forms 



Here B does not induce an inner product because the second term is not symmetric. Conse- 
quently the best approximation property in the energy norm is lost together with numerical 
stability, see [91 [33] . This has motivated the use of the Petrov-Galerkin method, where the 
combination of specially engineered test functions with added stabilization terms in the varia- 
tional formulation, has improved the numerical approximations considerably. The DPG method 
discussed here restores the best approximation property by the computation of a special set 
of test functions for given trial functions that produce a symmetric positive-definite algebraic 
system when substituted in the bilinear form. The method can be explored from many dif- 
ferent perspectives. It may be viewed as an automatic apparatus designed for guaranteeing 
the famous inf-sup condition of the Babuska-Brezzi theory of discrete stability, see [HE], or it 
may be interpreted as a generalized least squares method, see [T7]. We treat the methodology 
as a numerical construction which generates a true inner-product structure for any well-posed 
variational problem. 



B(w,u) = £(w), VweW 



(2) 




(3) 
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2.2 Petrov-Galerkin Method with Optimal Test Functions 

The goal is to find a linear operator T which maps the trial function space U to the test function 
space W such that 

<B(Tv, u) = A{v, u) Vv.ueW (4) 

where A(v, u) is a symmetric coercive bilinear form on U x U. If {ei, e2, . . . , e„} denotes a set 
of trial functions that spans a subspace U n of U, then the use of the test functions Tej in the 
Petrov-Galerkin method 

B(Te h u n ) = £(Tei), i = l,...,n (5) 

defines u n as the orthogonal projection of the exact solution u onto the trial space U n with 
respect to the inner product A(-, ■): 

.A(v n , u - u„) = Vv n G U n 

In other words, u n is the best approximation of u in U n in the corresponding 'energy' norm, 
that is, 

|||u-U n ||| W < |||U-V||| W Vv G Li n 

where 

Ill u |||w = V A(u, u) Vu G hi 

A canonical way to define the trial-to-test operator T : U — > W is to define Tu G W through 
the variational equality 

(w,Tu) w = S(w,u) Vw G W (6) 

so that we have 

B(Tv,u) = (Tv,Tu) w 

and thus achieve our goal with 



^(v,u) = (Tv,Tu) w , |||u||| w = ||Tu||w (7) 



2.3 Localization 



For each trial function ej, i = 1, . . . , n, Equation ^ constitutes an auxiliary, infinite-dimensional 
variational problem associated with the test function space W. The standard Galerkin method 
may be used to approximate solutions to Equation (jGb- However, for typical variational formula- 
tions, like for the one corresponding to the forms of ([3j), Equation ^ is from the computational 
perspective at least as 'heavy' as the original problem. 

The new insight of Demkowicz and Gopalakrishnan, see [16] , was to formulate the variational 
problem ^ in an ultra-weak hybrid form in such a way that the test function space W becomes 
fully discontinuous on a partitioning of Q into elements {Ki, K 2 , . . . , K^}: 

W = {w : w|^G W(K) for every element K} 

As long as the inner product (-, -) w does not couple test functions defined on individual elements, 
the auxiliary problems defined by Equation ^ can be solved locally at the element level. The 
localized version of the problem can be written in terms of the broken forms 

(v)w = E(v)w(jo. B(v) = £ B *(v) 

K K 
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which allows the global problem ^ to be written as finding Tu G W(K) such that 

(w, Tu) w(K) = B K {w, u) Vw G W{K) (8) 

This problem may then be solved approximately using the standard Galerkin method in an 
enriched finite element space Wh(K) C W(K), that is, find Tu G Wn(K) such that 

(w,fu) ww = S Jf (w,u) VwGW ft (if) (9) 
2.4 Well-posedness and Robustness 

The above construction can be carried out provided that the original variational problem ^ 
is well-posed, that is, provided that the bilinear form B(w,u) satisfies the conditions of the 
Babuska-Necas-Nirenberg Theorem, which are (see [TJ and [6| IT3]) 

B(w,u) < Ci||w||w||u|iw Vw g W, ueU (10) 

sup ^ W ' u) > C 2 ||w||w (11) 
ueu \\ u \\u 

sup ^> > ftllulk, (12) 
wew ll w ll>v 

The first condition implies that the right hand side of Equation (|6| is a bounded linear functional 
in W for every u6W. Therefore Tu G W exists and is unique. 

The second condition guarantees that the mapping T is surjective, that is, every w G W 
can be expressed in the form w = Tu for some ugW. Namely, if this was not the case, there 
would exist a non-zero w G W such that (w, Tu)w = for all u6W. But then also 

£(w,u) = VuGW 

by Equation ^ so that Equation (11) implies w = 0, which is a contradiction. 



Finally, the coercivity of the bilinear form A follows from condition (12). To see this, we 



recall that the norm of an element in a Hilbert space H can be computed by duality as 

ii ii ( v > u )« / 1Q \ 

\\u\\ n = sup — — j-; — (13) 

so that the energy norm can be expressed as 

(w, Tu) w B(w,u) 

|||u|||^ = sup — j-j — r: = sup 



?ew ll w llw wew Il w llw 



Now Equation (12) implies that 

^(u,u) > C 3 2 ||u||2 VuGW 

The characteristics of the energy norm depend on the values of the constants C\ and C3 
which in turn typically depend on the physical parameters and the problem geometry. As a 
result the Petrov- Galerkin method ^ may not be robust, that is, uniformly accurate with 
respect to the parameter values of a particular physical problem that is being modeled, see [2] . 
However, the parametric dependence can be removed completely from the bilinear form A(-, •) 
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by using a special norm for the test function space in (|6j). This norm is referred to as the 
optimal test space norm in the DPG literature, see for instance [48J, and can be defined as 



w w 



|T*w| 



u 



(14) 



where T* : W — > U is a linear mapping such that 



(T*w,u) w = B(w,u) VueW 



Equation (14) defines a proper norm in W again under the assumption that the variational 



problem is well-posed. Namely, the utilization of the duality argument (13), and conditions 



(10 )-(H) show that 



C 2 ||w||>v < ||T*w|L < Ci||w| 



w 



that is, 



vv 



is an equivalent norm to 



w- 



duality argument (13) yields 



Finally, condition ( 12 ) implies that T* is surjective so that a subsequent application of the 



l u lllw 



sup 



B(w,u) 

|w|||w 



sup 



(T*w,u 



u 



|T*w| 



u 



sup 



ll v llw 



\ u \\u 



(15) 



Thus, the introduction of (14) allows us to obtain a problem for each trial function which 



yields a robust method (convergent in the \ \ ■ \ \u norm as stated in (15)). This guaranteed best 
approximation property is what leads us to call the norm induced by Equation (14) optimal. 



As we will observe in the proceeding Section, the ultra-weak approach allows the expression of 
the optimal test space norm to be deduced directly from the bilinear form. Additionally, the 
leading terms can be expressed in a closed, localizable form. 



3 DPG Method for the Convection-Diffusion Problem 

In this section we follow [161 EE] and derive an explicit DPG variational formulation of the 
convection-diffusion equation. The ultra-weak DPG variational form of the problem is associ- 
ated initially to a partitioning of Q into elements denoted by flh- We will refer to the collection 
of all edges in the mesh as dQh = Uk9K. We will use the standard notation where L 2 (Q) 
and L 2 (Q) denote the spaces of square-integrable scalar- and vector- valued functions over Q, 
respectively. Moreover, iY 1 (f2) and H(div,Q) stand for the subspaces consisting of functions 
with square-integrable derivatives: 

H\n) = {ve L 2 (n) : Vv G L 2 (n)} 

h (div, n) = {r eL 2 (n) : v • r g L 2 (n)} 

3.1 Derivation of the Hybrid Ultra- Weak Variational Form 

The procedure begins by rewriting the second-order partial differential equation in as a pair 
of first-order equations 

e _1 cr = -Vtt , 

16 

V ■ cr + a ■ Vm = / 

Here we have chosen to define cr as the diffusive flux in order to be consistent with the above 
cited references. Other option would be to include the advective flux aw in the first equation 
defining cr as in our earlier works (121 H3]. In our experience, the difference between the two 
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approaches appears to be mostly cosmetic. In this work we will assume that the ambient fluid 
is incompressible so that V • a = and that the problem is formulated in a dimensionless form 
scaled such that lal ~ 1. 



Upon multiplying the first and second equations in (16) by a vector- valued test function 
r and a scalar-valued test function v, respectively, and integrating by parts element-wise, 
we obtain the following bilinear and linear forms corresponding to the abstract variational 
statement ^ 

B(w, u) = (e~V - Vv, cr)n h - (V • r + a • Vt> , u)n h + (t ■ n, u) 9 n h + (a n , v} dQh 

C(w) = (v,f) nh [ } 

Here n denotes the outward unit normal on each element edge and a n corresponds to the (outer) 
normal component of the total flux (cr + an) ■ n. Moreover, the test function is denoted by 
w = (r,f) and the trial function by u = (cr, u, u, a n ). For smooth functions, the notations 
(-, -)n h and (-, -)an h stand for the element-wise computed inner products 

(f,9h h = £ / fgdQ, (f,g) d n h = £ / fd ds 



The element integrals (-,-)n h in Equation (17) make sense provided that ct,u,t,v as well 



as V • r and Vt> are square integrable over each K. Thus, the test function space is defined as 
the broken space 

w = H(dw, n h ) x H\n h ) (is) 

where 

Jf(div, Q h ) = {t e L 2 (n) : V • r e L 2 (K) VK e Q h } 

H 1 ^) = {u6 L 2 (Q) : Vv G L 2 (K) \/K G Q h } 

In the general case, the interface action (-, -)an h in ( fl7| ) must be interpreted as an appropriate 
duality pairing between suitable generalized function spaces defined on dK. The right setting 
can be deduced from the theory of trace operators naturally associated to the spaces H (div, K) 
and H l (K). These spaces are formally denoted by H~ l l 2 (dK) and H l l 2 (dK) for the normal 
trace r • n \q K and the standard trace v\qk, respectively 

This leads to the definition of the trial function space as 

u = L 2 (n) x L 2 (n) x H^ 2 (dn h ) x H- l ' 2 {dn h ) (19) 

where the spaces for the interface variables u and a n are defined as 

H 1 J 2 (dQ h ) = {v\ dnh : veH^Q)} 
H-V 2 (dCl h ) = {77 U h : 77 G H(div,Q)} 

Here Hp(Q) refers to the subspace of if 1 (fi) encompassing the Dirichlet boundary condition 
u = g on dQ. The topology of the numerical trace spaces can be characterized for instance by 
using the extensions with minimal norm: 



{inf ||f||,H-i(n) : v G H 1 ^) such that v \an h = u] 



\\o n \\H-V\8n h ) = ( inf ll»7llff(div,n) : V e JET(div,n) such that r) ■ n\ dUh = a n } 
More details regarding these definitions can be found from [18]. For our purposes it suffices 

1 /2 

to note that piecewise polynomial functions in if D (dflh) must be globally continuous on dVth 
whereas discontinuities at the vertices of VL h are allowed in the space H~ 1 / 2 {dVlh)- 
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3.2 Finite Element Trial Spaces 

We refer to Pt as the set of polynomials of degree lower or equal to t and to the corresponding 
tensor product family by Q Px , Py = P Px x P Py and set Q p = Q PtP . The trial space is defined then 

as 

U n = S 2 h x S h x M h x N h 

where 

Sh = {v E L 2 (tt) : v\k G Q p for every element K} 
M h = {v e Hp(Q) : v\e G P p+ i for every edge E} 
N h = {v : v\e G P p for every edge £"} 

In this work we use the Bernstein polynomials as a basis for P t . The t + 1 basis polynomials 
are defined as 

Bi, t {x) = Qar^l-x)*"*, i = 0, 1, . . . , * 

3.3 Localization in the Convection-Diffusion Problem 



Upon selection of a norm 1 1 • | = (-, -)yy for the broken test function space W defined in ( 18 ), a 
precise form of the localized auxiliary problem ^ is obtained. This problem is solved approxi- 
mately using the standard Galerkin method and an enriched finite element space defined on each 
element K. We will associate this space to a sub-partitioning of K into elements {K±, . . . , K^}- 
More precisely, we define an ff(div, K) x H 1 (i^)-conforming finite element space for the 
test functions as 

Wn(K) = r h xv- h 

where 

T~ h = {t G if (div, K) : r|^G Qp+i,p x Qp^+i for every sub-element K} 

Vj t = {v G H l (K) : f|^-G for every sub-element K} 

Here we denote p = p + Ap so that the level of enrichment is controlled by both the number of 
sub-elements iV and the increment in the order of approximation Ap. In what follows, we will 
describe three different choices of the test space norm and discuss appropriate enrichments for 
resolving the corresponding test functions. 

The Standard Test Space Norm 

A norm known as the standard norm, or mathematician's norm, may be constructed by com- 
posing the natural norms of the spaces H(div,Qh) and H l (Qh) as in 

l|w||! N = ||V • r\\l h + ||r||^ + \\Vv\\l h + \\v\\l h 

This choice leads to two decoupled, standard variational problems for r and v over each element 
K which can be solved sufficiently accurately by using an enriched finite element space with 
Ap = 2 and N = 1, see [IS]. 

The Weighted Test Space Norm 

While the standard norm provides a local problem that is easy to solve, it leads to approxi- 
mations of u that underestimate the true values considerably. It was determined in [19], that 
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the use of a special weighting of the norm near the inflow boundary could improve the accu- 
racy. The original version of the weighted norm, as introduced in [19], is based directly on the 
standard norm and reads 

II w IIwn = IIV • r\\l hfi + \\r\\l h ^ + \\Vv\\l htP + |M|^3 



where /3 is a piecewise constant norm weight such that 

0(x) = 



7, if d(x,r~) < 5 & d(x,r+) > 5 
1, otherwise 



where d(x, T ) and d(x, T + ) stand for the (normal) distances between the point x and the 
inflow and outflow/no-flow boundaries, respectively. These are defined as 

r~ = {x e dQ : a • n < 0} 
T+ = {x G dfl : a • n > 0} 

It should be noted that different refined variants of the weighted norm have been recently 
introduced and studied in |20j . These refinements have resemblance to the quasi-optimal norm 
introduced next. 



The Quasi- Optimal Test Space Norm 

While the weighted test space norm resolves some of the deficiencies of the standard norm, the 
approach is not optimal in the sense of the DPG theory as outlined in Section [274} In contrast 
to the weighted test space norm, we describe here the optimal test space norm ( |14| ). However, 
the true optimal norm leads to a series of global problems so that some numerical modifications 
are necessary to make the approach suitable for practical computations. This way we arrive to 
a quasi-optimal test space norm which retains the essential features of the optimal test space 
norm but is localizable. 



Let the trial space norm be the natural norm arising from the definition (19) 



>. 1/2 

1 1 2 II'"' 1 1 ^ 1 1 2 11^112 I 

l u llw = i IMlLa(n) + IMIi 2 (n) + IMIizV2(anfc) + \\ a n\\H-v 2 (dn h )j 



Thus, using the definition (17) of the bilinear form, we can express the optimal test space norm 



(14) as 



HWw = He l r-Vv\\l + ||V-T + a. Vv\\l + llfr •nl||| 0t + ||M 



h 



where 



[T-n]||an h = sup — — , ||LvnJ||an fc = sup 



The last two terms can be made more precise, see |18j . but the important point is that they 
involve the 'jumps' of r and v along the inter-element boundaries and therefore prevent the 
element-wise computation of the optimal test functions. However, a localizable variant of the 
optimal norm may obtained simply by dropping these terms and adding an integral contribution 
of v, so as to remove the 'zero-energy mode' (t,v) = (0, const.). This leads us to the so-called 
quasi- optimal test space norm 



w 



Iqon = Ik ^-Vvll^ + HV-T + a-Vvll^+aillTll^ + aalHI^ (20) 
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where oci > and «2 > are numerical regularization parameters to be chosen. Notice that 
while the addition of a contribution of ||t||q is not necessary, we have observed that a proper 
regularization of r improves the accuracy of the final approximation. In our algorithm we 
choose a% = e~ 3 / 2 and a 2 = 1. 



A critical issue in the application of the quasi-optimal norm (20) is that the auxiliary 
problems for the test functions become singularly perturbed. This was explored in the one- 
dimensional case in [121H3] where a specially designed sub-grid mesh was employed to resolve the 
potential boundary layers in the test functions. To extend the idea to two spatial dimensions, 
we observe that the strong form of the problem (|8| associated to the quasi-optimal test space 



norm (|20|) takes the form 

-V(V ■ r) + (e~ 2 + a x )r - V(a • Vv) - e^Vv = e~ x a + Vm 



^ _ ^ . _ ^ (21) 



-Av - a ■ V(a • Vv) + a 2 v - a ■ V(V ■ r) + e X V ■ r = V • cr + a ■ Vu 



in every K. Notice that the system (21) corresponds to the Euler-Lagrange equations for the 
variational problem ^ and is accompanied with inhomogeneous natural boundary conditions 
involving also the interface variables a n and u. 

While a detailed regularity and asymptotic analysis of the system ( |2~T| is out of the scope of 
this paper, a preliminary examination of the system when K = (0, 1) x (0, 1) and a = (1, 0) or 
a = (0, 1) reveals that the solution may exhibit regular, exponential boundary layers behaving 
like 



e -*a>i e^ x ^ Xl \ e~ x ' X2 , e - x< ^- x ^) 



where the characteristic exponent is given by A = e _1 \/l + a\e 2 . Notice that when a\ = e 3//2 , 
we have A = e _1 + 0(e' 1 ^ 2 ). To account for these features in the automatic computation of the 
optimal test functions, we employ a layer-adapted Shishkin mesh on each K as shown in Figure 
[T] Such meshes are commonly used in the approximation of singularly perturbed problems, see 
[37J H6]. We follow a strategy which consists of adding one layer of 'needle' elements of width 
0(pe) near the boundary and has been proven to be optimal in the context of the hp- version 
of FEM for react ion- diffusion equations, see [101 H5| Wf\ . 

We emphasize that the sub-grid discretization is needed only when using the quasi-optimal 
test space norm. Failure to address the fine-scale features of the corresponding test functions 
(here using the sub-grid) degrades the accuracy of the final solution and dissipates all benefits 
of the quasi-optimal test norm. 



4 Algorithmic Considerations 

At the high level, the DPG method with optimal test functions is similar to the standard finite 
element method detailed in Algorithm [Tj Initially, the mesh information is input and mappings 
from the local stiffness matrix to the global matrix are determined. After this, we loop over the 
elements and form the element contributions to the global stiffness matrix and load vector. The 
main conceptual difference between the DPG method and standard FEM is the following: the 
local stiffness matrix will be computed as a by-product of computing the optimal test functions. 
Another difference is that, because of the hybridization, the trial spaces are non-standard and 
require some specialized handling. The remaining portion of this section will describe both 
differences in more detail. For a review of standard FEM implementation, see [TU 130] • 
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Figure 1: A 3 x 3 Shishkin mesh to resolve the boundary layers of the optimal test functions. 



Algorithm 1 The DPG finite element method with optimal test functions 

1: Initialize mesh 

2: Compute all local to global maps 

3: for all elements in mesh do 

4: Compute local contributions K^F^ using OptimalTestFunctions (Algorithm [2]) 

5: Assemble K<-^K and Fg — > F 

6: Apply boundary conditions 

7: end for 

8: Solve Ku = F 
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4.1 Selection of the Trial Spaces 

The space of trial functions in the DPG method consists of basis functions used to discretize 
variables on the element interiors (Sh in our application), as well on the edges (M h , N h ). Notice 
that on edges we must span spaces where corner degrees of freedom are assembled with those 
of neighboring edges (Mh) as well as spaces where corner degrees of freedom are independent 

What this means from the implementation standpoint, is that degrees of freedom from 
each space must be given unique numbers which correspond to their location in the global 
stiffness matrix. This is not straightforward because the basis functions are defined over varying 
topological objects (element interiors/edges) and discretize different variables, a consequence of 
using a first-order system. A flexible finite element framework is needed to address these issues. 
It is the authors' experience that even in the context of rectangular domains and structured 
grid meshes that this is a non-trivial task which is not immediately appreciated when digesting 
the DPG framework. 

4.2 Computation of the Optimal Test Functions 

The optimal test function problem ^ is solved locally at the element level. By specifying a trial 
function e U n one can compute the corresponding optimal test function w = Te^ 6 Wn(K) 
using the standard Galerkin method to solve 

(8w, w) w(J0 = B K (8w, e 4 ) V5w e W n (K) c W(K) 

Notice that since the trial function is given, the bilinear form i3^(5w,ej) induces a linear 
functional in W(K) and is used to form the right-hand side of the optimal test function problem. 
Notice that each trial function will lead to a different right-hand side in the above problem. 
We will denote by F opt the matrix whose columns are the load vectors corresponding to each 
trial function. The bilinear form on the left hand side is defined by the selected test space inner 
product and is independent of the trial functions. We will denote by K opt the corresponding 
stiffness matrix. The columns of the solution U opt to the system K opt U opt = F opt represent the 
approximated optimal test functions. 

Algorithm [2] describes the procedure used to form and solve this linear system. The al- 
gorithm is a function of the trial space U n , the enriched finite element space for the optimal 
test functions Wa(K), the chosen optimal test function inner product (-,-)w(k), and the bi- 
linear and linear forms Bk{-, •) and £(•) arising from the problem setup. At this point, we 
remind the reader that the optimal test functions are vector-valued. When forming the matrix 
K op t, the choice of optimal test space norm will produce a block matrix structure depending 
on which components of optimal test functions interact with each other. Recall that in our 
application we have three components 5w = (8ti,8t 2 ,Sv) and w = (r 1 ,r 2 ,f). We will denote 
the corresponding blocks of the matrix K opt by [K opt ]y, i,j = 1, 2, 3. 

The particular form of the matrix depends on the manner in which the inner product is 
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defined. If we choose the standard norm (SN), the blocks are defined as 



[Koptln = (<5ti,i,Ti,i)a- + {St^t^k 
[K opt ] 12 = (5T hl ,r 2 , 2 )K 



[K pt] 13 = 



[K op t] 21 = (5r 2)2 , Ti,i)a- 

[K opt ] 22 = (St 2;2 ,t 2;2 ) k + (5t 2 ,t 2 ) k 

[K op t] 23 = 

[Ko P t] 31 = 
[K op t] 32 = 

[K pt] 3 3 = (5v :1 ,v tl ) K + (5v i2 ,v i2 ) K + (5v,v) K 

where (•, -)k denotes the L 2 inner product computed over an element K and subscripted comma 
indicates differentiation with respect to the coordinate indicated by the subscript. If the 
weighted norm (WN) is used, the structure of the matrix is identical but the L 2 inner product 
is replaced by the appropriately weighted inner product. However, if we choose to work with 
the quasi-optimal norm (QON), the matrix K opt becomes full: 

[Koptln = (^1,1 > r i,i)^ + ( e ~ 2 + «i)( 5r i> t i)k 

[K op t] 12 — (^"1,1) r 2,2)K 

[K opt ] 13 = -e^^Ti^Jifi (<fri,i,aif,i + a 2 v )2 ) K 
[K op t] 2 i = (^ r 2,2,n,i)x 

[K opt ] 22 = (<Jr 2 , 2 , T"2,2)x + (e -2 + ai)(^T 2 , r 2 ) x 
[K opt ] 23 = -€~ 1 (8t 2 ,v :2 ) k + (ST 2 ,2,aiv tl + a 2 v >2 ) K 
[K opt ] 31 = -e _1 (<Ju i, ri)^ + (aify i + 02^,2, 
[K opt ] 32 = -e _1 (<k> )2 , r 2 ) x + + a 2 8v :2 , T 2t2 ) K 

[K opt ] 33 = (1 + al)(5v A ,v A ) K + (aia 2 5v t i, v j2 )k 
+ (1 + aj)(5v :2 , v :2 ) K + (a 1 a 2 6v t 2 , 

The rows of the right hand side F opt follow the same block structure as K opt . The columns 
correspond to different trial function components in u = (<7i, a 2 , u, u, a n ). We will denote the 
corresponding blocks of F opt by [F opt ]jj, i — 1,2,3, j — 1, ... ,5. The actual form of the matrix 
is determined by the bilinear form B(-,-): 



F op t — 



(e _1 5ri,(Ti)^ -(Sr 1A ,u) K (ST 1 n 1 ,u) dK 

{e~ l ST 2 ,a 2 ) K -{5t 2 , 2 ,u) k (Sr 2 n 2 ,u) dK 

-(<Ju,i,<Ti)a- -(8v, 2 ,cr 2 ) K -(a 1 5v i i + a 2 Sv t2 ,u) K (Sv,a n ) dK 



where (-,-)dK denotes the standard L 2 inner product computed over the element boundary dK 
and n\,n 2 denote the components of the unit outward normal to dK. 

The blocks in the matrices K opt and F opt are formed in the standard finite element fashion. 
Since the matrix K opt is mostly dense, we suggest the use of a direct solver with support for 
multiple right-hand sides, such as LAPACK's DGESV. 

Standard wisdom would then suggest that one could compute the local contributions 
in the natural way, by first looping over the space of optimal test functions and subsequently 
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the corresponding trial functions as in the Petrov-Galerkin method. However, the local matrix 
contributions may be computed directly by a matrix-matrix product of = Uo pt F opt . This 
can been seen most easily by writing the approximate optimal test function as 

h 

Tej = y^[U op t]fctefc 

k=l 

where {e 1; . . . ,e„} stands for the basis of the enriched finite element space Wa(K). Then it 
follows that 

h h 

[K £ ]ij = B(Tei,ej) = ^[IJ opt \kiB(e k , e,-) = ^[U op t]jfci[F opt ] fej 

k=l k=l 

so that 

\K(\ij = [U^Fopt]^- 

The local contribution to the problem right-hand side F;, however, must be assembled in the 
traditional way. From the implementation point of view, it is advantageous to think of DPG 
as a routine which computes a special local matrix which then must be globally assembled. 



Algorithm 2 Computation of the optimal test functions 



function OptimalTestFunctions(W„, V\4(A), (•, -)w(k), B k (-, •)> £(')) 
for all components 5wi of 5w do 
for all components Wj of w do 

if Swi interacts with Wj as defined by (•, -)w{k) then 

Assemble [Kopt]^- block of K opt 
end if 
end for 

for all components Uj of u do 

if 5wi interacts with 11 j clS defined by Bk{-, •) then 

Assemble [F opt ]jj block of F opt 
end if 
end for 
end for 

Solve K opt U op t = F opt 
K * = Uo pt F opt 



Assemble F^ by £(■) 
Return K^, F^ 
end function 



5 Numerical Results 

In this section we present numerical results which compare the performance of each of the 
optimal test space norms. While the DPG solution consists of the field variables cr,u and 
the element interface variables <r n ,-u, our visualizations will show only the u component of the 
total solution. In our experience it is helpful to visualize all solution variables, particularly 
while debugging, but the extra plots do not contribute here to a better understanding of the 



14 



performance of the method. However, we will also show numerical convergence results in terms 
of the L2 norm for u and cr which includes contributions from all solution variables. 

In all of our model problems the computational domain is taken to be the unit square 
Q = (0, 1) x (0, 1). The problems feature different Dirichlet boundary conditions and different 
orientations of the advection vector a with respect to the mesh. For each problem, we compare 
the coarse-mesh accuracy of the scheme for different values of the diffusivity (e = 10~\ i = 
2, 4, 6) and different test space norm. We shall use the abbreviation SN for the Standard Norm, 



WN for the Weighted Norm and QON for the Quasi-Optimal Norm, see Section |3.3| to recall 
their definitions. 

In the test problems we set p — 1 corresponding to discontinuous, piecewise bilinear rep- 
resentation of the field variables cr and u on the domain Q, discontinuous piecewise linear 
representation of the numerical flux a n and continuous, piecewise quadratic representation of 
the numerical trace u on the 'skeleton' dflh- The trial spaces are constructed on meshes of ten 
by ten elements. The optimal test functions are approximated using a polynomial enrichment 
value of Ap = 2. 

The Eriksson-Johnson problem 

Eriksson and Johnson consider in [23] the equation 

U i — eAu = f in Q 

representing the convection-diffusion equation with the flow velocity a = (1,0). The problem 
can be approached analytically by separation of variables. Assuming homogeneous Dirichlet 
boundary conditions at x = 1, y = 0, and y = 1 and the inflow data u(0,y) = sm(iiy), the 
solution becomes 

(l-s)x\ f(l + s)x-2s 



6XP I 26 ) ~ GXP V 2e 7 r- 

u(x,y) = — sin(7n/), s = y/T+ 47r 2 e 2 

1 — exp y — J 

and features a steep boundary layer of width e near the outflow boundary at x — 1, see Figure [2] 
for a schematics of the problem. As in all problem schematics, the dashed line will represent 
the location of the boundary layer and the shaded region represents the area weighted when 
using the weighted test space norm. 

Table [3] organizes a comparison of the norms for problems of increasing difficulty in terms 
of e versus the element size. A reference solution is included as a final row which in this case is 
the analytic solution. From the first row, one can observe the deficiency of the standard norm 
as applied to this problem. Significant undershooting of the exact solution is observed. This is 
important even in the context of adaptivity since error in the coarse scale solutions will lead to 
unwanted refinement. The success of the weighted norm in the second row also becomes clear. 
The rough character of the true solution is obtained which would drive adaptivity in the proper 
locations. However, the quasi-optimal norm is clearly superior in this case, capturing the fine 
scale features on the coarse mesh. 

Advection Skew to the Mesh: Continuous Dirichlet Data 

This problem has been studied in p^6l [19] and consists of solving Equation Q with the boundary 
conditions 

u(0,y) = l-y, u(x,0) = l-x, u(x, 1) = u(l,y) = 
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u = sm(iry) 




u = 

boundary layer 



weighted region 



u = 



u = 



Figure 2: Schematics of the Eriksson- Johnson problem. Approximate location of the boundary 
layer shown as a dashed line and the weighted region is shaded. 



and the advection vector a = (cos 6, sin 8), 9 £ [0°,90°]. The solution features a boundary 
layer of width e near the outflow boundary. We include a schematic of the problem in Figure [4] 
where the boundary layer is indicated by a dashed line and the weighted regions for the WN 
are shaded. For this problem, the reference solutions were computed using standard Galerkin 
FEM on a fine mesh which fully resolves the boundary layer. We were not able to compute 
such a solution for the problem with e = 10~ 6 although one can compare with the solution for 
e = 10~ 4 as the difference in this case is visually indistinguishable. 

We observe similar performance of the SN and WN as compared to the previous problem 
@. The QON in this case captures the main features of the solution well. However, while the 
WN exhibits a local overshoot in the solution for the e = 10~ 6 case, the solution value for the 
QON is underestimated across the domain. This occurs probably because the test functions 
are not sufficiently well resolved by our sub-grid discretization, see comments ahead. 



Advection Skew to the Mesh: Discontinuous Dirichlet Data 

Our last problem is the standard benchmark test for the convection-diffusion equation in the 
engineering literature. The advection vector is the same as in the previous problem, but the 
inflow data is taken to be piecewise constant: 



u(x,0) = l, u(y,0) 



1, < y < 0.2 
0.2 < y < 1 



The problem schematic is shown in Figure [6] and the results in Figure [7j We emphasize here 
a feature of QON: if the mesh does not resolve the boundary layer, its effect is included but 
the layer is not resolved. This is in contrast to the SN and WN, where we always observe a 
boundary layer the width of one element. As in the previous problem, the reference solution 
has been obtained on a fine mesh using the standard Galerkin FEM. 



5.1 Norm Comparison: Mesh Refinement 

We study the convergence of the different DPG schemes in the advection skew to the mesh 
problem with continuous Dirichlet data under uniform h and p refinements of the trial space. 
We compare the error of u and cr with respect to the reference solution measured in the L 2 norm 
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Figure 3: Color map plots of u in the Eriksson- Johnson problem with an analytic solution. 
Color range is uniformly set to [0, 1] and the minimum and maximum values of the solution are 
shown underneath each plot. 

e = KT 2 e = 10" 4 e = 10" 6 
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u = 



u = l-y 



■ 


boundary layer , 




/ i 




a 1 




weighted region i 




/ 1 


/ 1 


u — 1 — X 



Figure 4: Schematics of the advection skew to the mesh problem with continuous Dirichlet 
data. Approximate location of the boundary layer shown as a dashed line and the weighted 
region is shaded. 

when e = 10~ 2 . The purpose of this study is to show quantitatively the improved accuracy of 
the QON solution. 

Results from this study are presented in Tables [T] and |2j Note that in each table, the 
reference solution comes from the standard Galerkin FEM with a fine mesh of rectangular 
bilinear elements. The values of er ref were obtained by a linear interpolation of the coefficients 
determined from taking finite differences of the solution u Te {. In both refinement strategies and 
for both u and cr, the QON provides better approximation to the reference solution, particularly 
when the mesh is coarse. 



5.2 The Resolution of the Test Functions Corresponding to the 
Quasi-optimal Test Space Norm 

The best approximation property of the DPG method hinges on the sufficient resolution of the 
optimal test functions. In the context of the quasi-optimal norm this is non-trivial because of the 



singular perturbation character of the auxiliary problem (21). While the subgrid discretization 
we have proposed is an effort to address this issue, the link between the accuracy of the test 
functions and the final solution is still unclear. This is partially due to the fact that the 
auxiliary problem is, to the best of our knowledge, unknown in the literature. Therefore, we 
make a numerical convergence assessment of the optimal test functions. 

We restrict ourselves to p-refinements on the subgrid and measure the approximate relative 
error of the test functions in the quasi-optimal test space norm ||| • |||w(-ftr) as m 



wP +1 - w 



P \\\w(K) 



W 



p+1 



\\\w(K) 



This we carry out on a single square element K of the ten by ten mesh used earlier in the 
advection skew to the mesh problems. Because the norm ||| • is the usual energy norm 

for the auxiliary problems, its computation reduces to linear algebra. Recall that the columns 
of the matrix U op t represent the optimal test functions. Therefore the energy norm squared of 
the optimal test functions Wj is given by 

lll w ;|llvv(/0 = U opt (:,i) T K opt U opt = U opt (:, i) T F opt (:, i) = K e (i,i) 
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Figure 5: Color map plots of u in the advection skew to the mesh problem with continuous 
Dirichlet data. Color range is uniformly set to [0, 1] and the minimum and maximum values of 
the solution are shown underneath each plot. 

e = HT 2 e = 10" 4 e = 10" 6 
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Figure 6: Schematics of the advection skew to the mesh problem with discontinuous Dirichlet 
data. Approximate location of the boundary layer shown as a dashed line and the weighted 
region is shaded. 



Table 1: Convergence of u and cr under uniform h- and p-refinement in the advection skew to 
the mesh problem with continuous Dirichlet data at e = 1CT 2 . 



Urel — U n \\L 2 (n) \ l^rcf ~ & n\ \L 2 (Si) 



N 


SN 




WN 




QON 




SN 




WN 




QON 




5 


2.55x10" 


-i 


1.17x10" 


-i 


7.69x10" 


-•2 


4.28x10" 


-2 


4.15x10" 


-2 


3.38x10" 


-2 


10 


1.76x10" 


-i 


6.54x10" 


-2 


4.84x10" 


-2 


3.80x10" 


-2 


3.68x10" 


-2 


2.45x10" 


-2 


25 


6.26x10" 


-2 


2.20x10" 


-2 


1.43x10" 


-2 


2.37x10" 


-2 


2.33x10" 


-2 


9.56x10" 


-3 


50 


1.64x10" 


-2 


7.35x10" 


-3 


4.33x10" 


-3 


1.09x10" 


-2 


1.08x10" 


-2 


3.13x10" 


-3 


100 


3.24x10" 


-3 


2.18x10" 


-3 


1.17x10" 


-3 


3.76x10" 


-3 


3.74x10" 


-3 


8.59x10" 


-4 



Mref — Un\\L 2 (n) \ |°"ref ~ & n\ \L 2 (Si) 



V 


SN 




WN 




QON 




SN 




WN 




QON 




1 


2.55x10" 


-i 


1.17x10" 


-i 


7.69x10" 


-2 


4.28x10" 


-2 


4.15x10" 


-2 


3.38x10" 


-2 


2 


1.49x10" 


-i 


5.56x10" 


-2 


3.97x10" 


-2 


3.45x10" 


-2 


3.29x10" 


-2 


2.36x10" 


-2 


3 


7.00x10" 


-2 


2.72x10" 


-2 


2.11x10" 


-2 


2.37x10" 


-2 


2.28x10" 


-2 


1.43x10" 


-2 


4 


2.58x10" 


-2 


1.30x10" 


-2 


1.07x10" 


-2 


1.40x10" 


-2 


1.37x10" 


-2 


7.66x10" 


-3 
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Figure 7: Color map plots of u in the advection skew to the mesh problem with discontinuous 
Dirichlet data. Color range is uniformly set to [0, 1] and the minimum and maximum values of 
the solution are shown underneath each plot. 

e = KT 2 e = 10" 4 e = 10" 6 
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Table 2: Convergence of u and cr under uniform h- and p-refinement in the advection skew to 
the mesh problem with continuous Dirichlet data at e = 10 -4 . 



Mref — U n \ \L 2 (n) \ 1 CT ref ~ C ra 1 1 L2 (Q) 



AT 


SN 




WN 




QON 




SN 




WN 




QON 




5 


3.73x10- 


-i 


1.72x10- 


-i 


7.05x10- 


-2 


5.23x10- 


-3 


5.23x10- 


-3 


5.22x10- 


-3 


10 


3.60x10- 


-i 


1.25x10- 


-i 


6.50x10- 


-2 


5.23x10- 


-3 


5.22x10- 


-3 


5.20x10" 


-3 


25 


3.51x10- 


-i 


9.12x10- 


-2 


6.61x10- 


-2 


5.22x10- 


-3 


5.21x10- 


-3 


5.16x10' 


-3 


50 


3.45x10- 


-i 


7.71x10- 


-2 


6.46x10- 


-2 


5.21x10- 


-3 


5.18x10- 


-3 


5.08x10- 


-3 


100 


3.36x10- 


-i 


6.70x10- 


-2 


6.15x10" 


-2 


5.17x10- 


-3 


5.13x10- 


-3 


4.91x10- 


-3 



| Pref — U n \ |L 2 (n) 

p SN WN QON 

1 3.73XHT 1 1.72X10" 1 YMxW 2 

2 3.56X10- 1 1.24X10- 1 3.44X10" 2 

3 3.52X10" 1 LOOxlO" 1 2.58xl0~ 2 

4 3.49X10" 1 9.16xl0~ 2 2.08xl0" 2 



1 1 C rcf — 0"n\\L 2 (d) 

SN WN QON 

5.23xl0~ 3 5.23xl0- 3 5.22xl0 -3 

5.23X10" 3 5.22xl0- 3 5.20X10" 3 

5.22X10" 3 5.20xl0- 3 5.17X10" 3 

5.20xl0~ 3 5.18xl0" 3 5.13xl0- 3 



where we have adopted the MATLAB notation for indexing a whole row or column of a matrix. 
In other words, the energy norm of the optimal test functions may be read directly from the 
diagonal of K^. 

Notice that when p = 1, there are 28 trial functions and optimal test functions. To simplify 
the representation of the results, we separate the optimal test functions into four groups asso- 
ciated to the trial functions for cr, u, u, and a n and compute the average of the errors for each 
group. We present the results in Tables [3] and [4] for e = 10~ 2 and e = 10 -4 where we have used 
e to indicate an average error and subscripts to refer to each group. We have also included the 
L 2 error with respect to the reference solution of u in the table. The results show that while 
the optimal test functions corresponding to the field variables cr and u are seemingly easy to 
approximate, those corresponding to the interface variables a n and u are problematic. As a 
matter of fact, the relative error level of the optimal test functions corresponding to the trace 
u barely decreases below the discretization error level for the field varible u in the trial space 
as Ap increases. This slow convergence explains the lull in the /i-convergence in Table [2] as well 
as the degradation of the global accuracy when e = 10 -6 in Figures [5] and [7j These phenomena 
are manifestations of the current practical limitations in the utilization of the quasi-optimal 
norm. 
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Table 3: Convergence of the optimal test functions in the advection skew to the mesh problems 
when e = 10~ 2 , N = 10 and p — 1. 



average relative approximate errors 

6/x ^7/, ^"IL &rr~ 



Ap 

1 2.16xl(r 3 4.85xl(T 3 



2 3.28xHT 4 

3 3.08xl0- 5 

4 2.52xl0- 6 

5 4.22X10" 7 

6 1.08xl(r 7 

7 3.44xl0- 8 

8 1.27xl(r 8 

9 5.22X10" 9 



7.38xl0- 4 
6.74xl0- 5 
5.34xl(T 6 
8.70xl0- 7 
2.20X10" 7 
6.92xl0- 8 
2.53xl0- 8 
1.04xl(r 8 



9.42X10" 2 
1.42X10" 2 
1.77xl0- 3 
2.15xl0- 4 
4.23xl(T 5 
1.12xl0- 5 
3.59xl0- 6 
1.32xl0- 6 
5.45xlQ- 7 



^ n 

1.13xl0- 3 
1.67xl0- 4 
1.57xl0- 5 
1.30xl0- 6 
2.18xl0- 7 
5.58xl0- 8 
1.77xl0- 8 
6.54xl0- 9 
2.69xl0- 9 



l^ref — W n ||x, 2 (n) 

4.62xl0- 2 
4.79xl0- 2 
4.84xl0- 2 
4.85xl0- 2 
4.85xl0- 2 
4.85xl0- 2 
4.85xl0- 2 
4.85xl0- 2 
4.85xl0~ 2 



Table 4: Convergence of the optimal test functions in the advection skew to the mesh problems 
when e = 10~ 4 , iV = 10 and p = 1. 



average relative approximate errors 



Ap 


















1 1 ^ref 1 1 


L 2 (Q) 


1 


2.07x10" 


-4 


5.40x10" 


-3 


4.54x10" 


-i 


5.36x10" 


-3 


9.15x10 


-2 


2 


1.51x10" 


-5 


3.00x10" 


-5 


3.25x10" 


-l 


7.65x10" 


-4 


7.22x10 


-2 


3 


8.17x10" 


-6 


1.68x10" 


-5 


1.48x10" 


-l 


3.99x10" 


-4 


6.50x10 


-2 


4 


3.10x10" 


-6 


6.22x10" 


-6 


4.67x10" 


-2 


1.50x10" 


-4 


6.32x10' 


-2 


5 


1.17x10" 


-6 


2.18x10" 


-6 


1.43x10" 


-2 


5.75x10" 


-5 


6.29x10 


-2 


6 


5.36x10- 


-7 


9.49x10" 


-7 


5.77x10" 


-3 


2.67x10" 


-5 


6.29x10" 


-2 


7 


2.97x10- 


-7 


5.02x10" 


-7 


3.56x10" 


-3 


1.49x10" 


-5 


6.29x10" 


-2 


8 


1.85x10" 


-7 


3.09x10" 


-7 


2.76x10" 


-3 


9.24x10" 


-6 


6.30x10" 


-2 


9 


1.24x10" 


-7 


2.04x10" 


-7 


2.39x10" 


-3 


6.15x10" 


-6 


6.31x10 


-2 
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6 Concluding Remarks 



We have studied the robustness of the discontinuous Petrov-Galerkin method with optimal test 
functions when applied to convection-diffusion problems. In particular, we have introduced a 
DPG scheme based on the quasi-optimal test norm and compared it with schemes based on other 
test space norms. Our numerical tests reveal that the quasi-optimal norm leads to improved 
coarse mesh accuracy as compared with the conventional DPG test space norms provided that 
the associated test functions are sufficiently well resolved. We have addressed the resolution 
of the optimal test functions by studying their convergence numerically. We have been able to 
find a compromise between accuracy and stability of the final DPG solution by combining a 
carefully designed element sub grid with numerical regularization of the test functions. 

The DPG method with optimal test functions originates from the functional analytic foun- 
dation of the finite element method, and is perhaps still far from being adopted into engineering 
practice. We have made an attempt to popularize the methodology by contributing a holis- 
tic approach which begins in the abstract framework and ends in detailed explanation of the 
algorithm from the implementation viewpoint. 
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